Bulletin of Mathematical Biology

We present a multiscale model of HIV that couples within-host viral dynamics with population-level transmission to capture the interplay between individual infection and epidemic spread. The model is structured by treatment age, allowing viral load to influence both infectiousness and progression to AIDS. The model is fitted using both clinical data (viral load and target cell counts of treated individuals) and epidemiological data (HIV incidence, diagnoses, and AIDS classifications). We derive the basic reproduction number and establish threshold conditions for the existence and stability of disease-free and endemic equilibria. Structural identifiability is assessed via input-output equations, showing that the multiscale model is identifiable when the initial number of treated individuals is set to zero and the AIDS death rate is assumed known. Parameters are estimated sequentially: within-host parameters are obtained via nonlinear mixed-effects modeling of clinical data, followed by estimation of population-level parameters using CDC surveillance data. Numerical simulations are performed using a finite-difference scheme with Picard iteration, and practical identifiability is evaluated via Monte Carlo simulations across varying noise levels. Results indicate that current strategies are unlikely to meet the 2030 targets, while increasing diagnosis rates and reducing transmission from diagnosed individuals could significantly alter epidemic trajectories. These findings highlight the importance of multiscale modeling and identifiability in informing effective HIV intervention strategies.

We propose a minimal mathematical framework to describe the go-or-grow dynamics of tumor cells comprising two phenotypically distinct populations. One population is migratory and undergoes constant diffusion, while the other proliferates in an oxygen-dependent manner. The local oxygen concentration governs transitions between these phenotypes. We then ask whether the time scale of cyclic hypoxia, where oxygen is spatially homogeneous but fluctuates between normoxic and hypoxic levels over time, alters the dynamics of the reduced generalist model in the fast phenotype-switching regime compared with the full go-or-grow model. We observe that the reduced generalist phenotype model fails to capture the go-or-grow dynamics under rapid oxygen fluctuations. Our results predict that the tumor growth rate increases in the generalist-phenotype model under rapid oxygen fluctuations compared with tumors of the proposed go-or-grow model. We then explore whether oxygen consumption, which leads to spatial heterogeneity in oxygen concentration, plays a role in distinguishing the dynamics of the reduced model from the proposed go-or-grow model. We demonstrate that the tumor in the proposed go-or-grow model exhibits a higher growth rate than the tumor in the reduced generalist-phenotype model under a slow oxygen consumption rate.

Fractional differential equations are increasingly used to represent temporal memory in tumour immune and pharmacokinetic models, but changing the derivative order without a consistent time scale can make comparisons difficult to interpret. We study a dimensionless, commensurate Caputo system for tumour, immune, normal cell, and drug states. We establish positivity, boundedness, the tumour invasion threshold, and local stability conditions. We then use a fractional Adams Bashforth Moulton predictor corrector method with an explicit mesh convergence study. Two conclusions correct potentially misleading interpretations of the model. First, the equilibrium locations are independent of the fractional order: decreasing the order changes transient time scales but does not by itself alter the model's asymptotic states. Second, after a unit bolus the drug equation has the exact solution C({tau})=E(-K{tau}). Its maximum is therefore the initial concentration, not a delayed peak. For 0<<1, the algebraic tail makes the infinite-horizon area under the concentration curve diverge. In the illustrative coupled simulations, the apparent late tumour reduction at small is consequently driven by persistent modelled drug exposure and must not be read as evidence of superior therapy. All parameters are dimensionless and illustrative; no patient data are fitted. Source code and generated tables are included to make every numerical statement reproducible.

Tumors consist of heterogeneous phenotypic cells, such as normoxic cells, which are highly proliferative, and hypoxic cells, which are less proliferative. Their phenotypic switching depends on tumor microenvironmental factors, such as oxygen and nutrient concentrations supplied by local blood vessels. However, during ongoing angiogenesis, the process of sprouting new blood vessels at the tumor site from pre-existing blood vessels, and how this phenotypic switching affects and impacts tumor growth, remains poorly understood. In this article, we formulate a mathematical model to elucidate the crosstalk between vasculature and tumor cellular heterogeneity during tumor progression. The model results show a strong agreement with the experimental data. Our simulation results demonstrate that ongoing angiogenesis increases tumor growth rate. In addition, we observe that the influence of hypoxic cells on phenotypic switching from normoxic to hypoxic is more pronounced than their influence on the transition from hypoxic to normoxic. Furthermore, we perform a global sensitivity analysis using the Sobol's method to assess the importance of the model's parameters. It highlights that the volume at which blood vessels attain half-maximal rate has the maximum effect on the model.

Tumor progression is often modeled as a passive response to external therapy or immune pressure, but tumor populations may also exhibit population-level regulation of proliferation and apoptosis. We develop a continuous-time Markov decision framework in which a controlled birth--death process represents a tumor population modulating the balance between proliferation and susceptibility to apoptosis in the presence of extrinsic death pressure. We examine threshold and quadratic costs, an unbounded linear reward, and constrained linear and quadratic formulations to determine how objective structure shapes optimal policies and induced population drift. Threshold and quadratic penalties generate restoring dynamics, with transitions from growth to suppression and regions of near-neutral drift associated with regulated or near-dormant behavior. An unbounded linear reward instead produces sustained or near-neutral growth without a restoring regime. Under constraints, a linear reward expands the region of positive drift as capacity increases, whereas a quadratic reward generates restoring, logistic-like drift around an interior population scale. These results show that regulated tumor dynamics depend on how growth incentives, extrinsic death pressure, penalties, and constraints scale with population size.

Much of the literature and intuition associated with mathematical epidemiology is driven by deterministic models, which are a reasonable assumption when the population size is large. Stochastic models, especially individual based models, are however considered vital when dealing with small population sizes, especially at times of invasion or extinction. The overwhelming majority of these models (both deterministic and stochastic) assume that the underlying parameters are fixed (or follow a regular seasonal pattern). Here, we consider an analytic framework for dealing with randomly varying parameters through the use of stochastic differential equations - thereby capturing the action of external noisy processes such as weather. In particular, we focus on when the transmission rate, {beta}, varies as the solution to a Cox-Ingersoll-Ross Model, such that {beta} is gamma distributed with autocorrelation. We consider the impact of this parameter variation on a stochastic version of the Susceptible-Infected-Recovered model, and for this 'double-stochastic' model show through simulation and analytical results that exogenous noise increases the impact of stochasticity, potentially leading to more early extinctions, wider variations in the number of cases at equilibrium, but that early growth rate can be faster or slower depending on the precise parameters.

Moment-closure methods provide a widely used approach for approximating low-order statistics of stochastic biochemical reaction networks, but the resulting closed equations may exhibit physically meaningless dynamics. Here, we systematically investigate the stability of second-order moment-closure approximations for a bursty gene expression model with nonlinear protein degradation. We consider the Gaussian, log-normal, and negative-binomial closures. By analyzing the resulting dynamical systems, we show that the Gaussian closure fails the proposed stability criterion for all parameter values, whereas both the log-normal and negative-binomial closures possess a unique physically accessible fixed point and are stable if and only if kB(1 + 2B) [&ge;] 1, where B is the mean burst size and k is the dimensionless protein synthesis rate. Moreover, this condition is automatically satisfied when the steady-state mean protein copy number is at least one, suggesting that it is mild in biologically relevant regimes. Numerical comparisons over a broad parameter range further show that the negative-binomial closure is the most accurate, while the Gaussian closure is the least accurate. We further develop higher-order negative-binomial moment-closure schemes and demonstrate that increasing the closure order substantially improves predictive accuracy. These results highlight the importance of considering dynamical stability, in addition to numerical accuracy, when assessing moment-closure approximations.

Ecological interactions play a crucial role in driving speciation and adaptive radiation. According to the adaptive dynamics theory of evolution in one-dimensional trait spaces, such evolutionary diversifications can occur when evolutionary branching points (convergence stable and evolutionarily unstable singular points) exist in the trait spaces. However, these evolutionary diversifications may occur across multiple traits at once, as exemplified by various within-species ecotypic diversifications involving multiple traits or genes. This paper attempt to extend this branching point condition to evolution in multi-dimensional trait spaces. Two criteria, termed branching possibility and branching inevitability, are developed to evaluate the likelihood of evolutionary branching in multi-dimensional trait spaces. Branching possibility guarantees the existence of evolutionary pathways that lead to branching (where deterministic parts of these pathways are described with the canonical equation of adaptive dynamics theory). Branching inevitability guarantees a steady progression of evolutionary branching (where the progression degree is described with a time-increasing function related to the local Lyapunov function, the potential function in an evolutionary potential game, and Fisher's fundamental theorem). Numerically simulated evolution in multi-dimensional trait spaces show that points having branching possibility or inevitability can induce evolutionary branching not only under asexual reproduction but also under sexual reproduction.

Biological populations have been evolving in multi-dimensional trait spaces. Ecological interactions can induce the evolutionary diversification of a biological population into distinct multiple populations, which is called evolutionary branching in adaptive dynamics theory. If possible mutations are restricted to certain directions due to developmental, physiological, or physical constraints, the resulting evolutionary dynamics may be restricted to subspaces with fewer dimensions than the original trait spaces. In general, such constraints are also evolving slowly, due to rare and slight mutational deviations from the constraint subspaces. In this sense, constraints are not strict. However, the method for describing non-strict constraints that is globally consistent has yet to be developed. This study proposes a general and flexible method for constructing a non-strict constraint subspace from the corresponding strict constraint subspace. This method enables efficient analysis of adaptive evolution including evolutionary branching under non-strict constraints in the framework of adaptive dynamics theory. Applying this method to a resource competition model shows that the shape of the carrying capacity distribution can be treated as a non-strict constraint curve, and that the likelihood of evolutionary branching is affected not only by the shape of the constraint curve but also by its slow evolution.

Objective: Extracellular electric fields critically influence neural activity through both exogenous neuromodulation and endogenous ephaptic coupling. While conventional cable models efficiently simulate membrane dynamics, they fail to capture bidirectional, field-mediated interactions self-consistently, and fully coupled volumetric methods require computationally prohibitive 3D meshing. We present Cable-BEM, a hybrid wire-kernel bidomain boundary element method designed to resolve these limitations. Approach: By analytically integrating boundary integral kernels around cylindrical neuronal compartments, Cable-BEM fully couples intracellular, extracellular, and membrane dynamics while strictly retaining the highly efficient 1D degrees of freedom of traditional cable equations. The system is advanced using a semi-implicit Crank-Nicolson scheme. To overcome the dense nature of the resulting integral operators, we implement an Adaptive Cross Approximation (ACA) and Hierarchical Off-Diagonal Low-Rank (HODLR) compression scheme. Main result: The solver was rigorously validated against full-surface bidomain boundary element method (BEM) reference implementations, demonstrating tight agreement in activation thresholds (within 1.3% relative error) across diverse stimulation geometries. The ACA-HODLR compression scheme achieved substantial memory footprint reductions-by a factor of up to 4.6 for large 225-cell networks- without sacrificing numerical accuracy. Furthermore, we utilized the framework to resolve subtle, distance-dependent ephaptic interactions, successfully demonstrating the progressive phase synchronization of biophysically realistic, multi-compartment Purkinje cells. Significance: Cable-BEM provides a computationally scalable, mesh-free framework that establishes a powerful and practical foundation for investigating complex field-mediated phenomena in large-scale, multicellular neuronal networks.

The global potential of a chemical reaction network has many applications and is closely related to the stochastic detailed balance. However, many fundamental questions concerning stochastic detailed balance remain unresolved, such as whether it depends on the system volume and how to construct new systems that satisfy it. In this paper, we show that stochastic detailed balance may depend on the system volume. We therefore introduce four types of stochastic detailed balance according to their dependence on volume and rate constants, and systematically investigate the relationships among them. Our results distinguish detailed balance arising from particular choices of volume and parameters from that enforced by network structure, and identify conditions under which detailed balance at one volume extends to all volumes. We further obtain a class of networks satisfying stochastic detailed balance for every volume and every positive choice of rate constants, and construct new systems whose global potentials exhibit double-well structures.

In many biological, biomedical and industrial systems, particles are transported via fluids through confined networks, in which clogging can disrupt function. However, we lack a predictive theoretical framework that couples particle transport, suspension rheology and network flow resistance. Here, we develop a model and solution algorithm for particle suspension flow in networks based on vessel-level continuum modelling and particle distribution at nodes connecting vessels. We apply the model to study transport of dense particle suspensions in minimal and physiological biological networks. A key feature of our model is the coupling between particle volume fraction and particle flux: in line with the physics of dense suspensions, each network branch, or vessel, possesses a local carrying capacity for particle transport at an intermediate particle fraction between zero and the maximum packing fraction. If this flux capacity is reached, the vessel becomes flux-limited and particles can accumulate in upstream branches, causing them to enter a high-particle-fraction, high-resistance state that we refer to as 'clogged'. We show that these vessel flux limitations lead to network-level redistribution of particles, which can cause widespread clogging and emergent network-scale heterogeneity. By varying network topology, we find that in some regimes increasing network connectivity does not improve transport: paradoxically, additional pathways can promote clogging and reduce network-level particle flux, analogous to classic results in traffic flow networks. Our results provide a minimal mechanistic framework that links suspension physics, network topology, and transport failure in complex flow networks.

Buruli ulcer (BU) is a neglected tropical skin disease caused by infection with Mycobacterium ulcerans. The mode of its transmission into the skin has not been fully clarified, but it is believed to involve mechanical transfer, perhaps by insect bites. M. ulcerans secretes a toxin called mycolactone which causes necrosis of the skin and formation of large skin lesions. Understanding the development of these skin lesions over time is crucial to understand the pathogenesis of BU. We formulate a mathematical model to describe the space-time dynamics of M. ulcerans bacteria, mycolactone and skin cell densities as well as fibrin deposition. We investigate the relative contributions to ulcer formation of direct cytotoxicity of mycolactone and ischaemia following the loss of coagulation control. Our results are consistent with existing biological evidence that higher mycolactone concentrations accelerate cell death. Furthermore, the model indicates that the rate of wound expansion is primarily controlled by the diffusion and net birth rate of bacteria, which sets the overall scaling for the speed of wound expansion, along with two dimensionless constants that define the shape of the trajectory of wound edge position against time, namely the ratio of the mycolactone diffusion coefficient to the bacterial diffusion coefficient, and the ratio of the mycolactone degradation rate to the bacterial net birth rate. When mycolactone diffuses significantly faster than bacteria spread or when mycolactone degrades more slowly than bacteria reproduce, the initial development of the wound is faster. In contrast, the long-term wound expansion speed appears to be relatively unchanged across a wide range of the dimensionless parameters.

Pseudoprogression in response to immune checkpoint inhibitor (ICI) therapy can be a complex challenge during treatment response evaluation. Erroneous classification of progression may result in stopping effective treatment for patients, while delayed confirmation of progression can prolong an ineffective treatment course. We develop a mechanistic mathematical model of tumor-immune dynamics under ICI treatment in which observed tumor volume includes contributions from both tumor burden and immune infiltration. We show analytically how pseudoprogression can arise from increased immune infiltration and proliferation as well as continued tumor growth during the ``ramp-up" phase. The model also incorporates circulating tumor DNA (ctDNA) dynamics, which we demonstrate are able to closely fit existing longitudinal tumor and ctDNA data and extract realistic parameter ranges to generate representative in silico data. Using model-generated cohorts, we train a random forest (RF) classifier using longitudinal tumor and ctDNA measurements to distinguish pseudoprogression and true progression, and demonstrate that it improves response classification relative to the Response Evaluation Criteria in Solid Tumors (RECIST 1.1) and provides classification earlier than the immunotherapy-specific iRECIST, even in the presence of tumor measurement noise. This work predicts that paired longitudinal tumor-ctDNA measurements may enable earlier discrimination of true progression from pseudoprogression, and motivates collection of imaging and ctDNA data synchronously during ICI therapy, especially in the setting of suspected disease progression.

Spatial population structure and seasonality are both central to the spread of many infectious diseases of plants, animals and humans. While seasonal forcing in transmission often plays an important role in epidemiological models of a wide range of infectious disease, and we now have some theoretical understanding of the dynamical impacts of spatial structure, the combined effects of these two ubiquitous processes has not been examined in detail. Here, we develop a novel model to explore the combined influence of spatial structure and temporal variability on disease dynamics. Spatial structure is represented using a lattice-based approach with near-neighbour interactions, while temporal variability is included through regular, seasonal, variation of the transmission rate. We use bifurcation analysis of a pair approximation of the full spatial model to identify the parameter regimes associated with qualitatively distinct dynamical behaviours. The model exhibits a remarkably wide range of complex dynamics, including limit cycles, quasi-periodic cycles, multi-year cycles, chaotic dynamics and bistability between these different states. In particular, complex dynamics occur when reproduction is predominantly local, with the dynamics depending critically on the amplitude of the seasonal transmission rate. We show how high transmission rates, high birth rates and in particular low recovery rates are requirements for complex dynamics. We predict that SI-type disease interactions in plant pathogen systems will show complex dynamics even with relatively global transmission dynamics.

The study of tissue dynamics has been stimulated during the last decades thanks to the use of quantitative descriptions, with the development of several theoretical and computational frameworks, many of them revolving around the notion of reaction-diffusion systems, eventually with additional structure variables beyond time and space. The use of structure variables can accommodate phenotypic traits. In this work, we study a family of competition models, where a given population depends on a resource (e.g. oxygen) and several populations are competing for it. Our quantitative description incorporates phenotypic traits and heterogeneity at the level of cell cycle variations, which influence replication rates via oxygen consumption. This enables us to replicate the fitness of specific subpopulations to environmental conditions (e.g. oxygen shortage or external influences). Using numerical simulations, we show that such models display dynamical pattern formation in the form of coupled travelling wave profiles that expand or retreat at the same wave speed. The full theoretical analysis of such dynamics is quite involved; to circumvent this difficulty, we introduce a quasi-stationary approximation for the resource dynamics. We find that this approximation can reproduce the overall behaviour very accurately, with the additional benefit of allowing theoretical treatment of the reduced model. In this way, we provide estimates on the wave speed which are numerically shown to be robust across a wide range of macroscopic parameters of the full model. The wave speeds are thus found to depend strongly on the proliferation rate of the fittest population, resembling a winner-takes-all dynamics.

Immune checkpoint blockade (ICB) therapy has revolutionized cancer treatment, though it is still effective in only about 20-40% of patients. To increase the efficacy of ICB therapy, combinations of ICB antibodies such as anti-PD1 and anticoagulants (e.g., thrombin inhibitors) have been studied in preclinical mouse models and clinical trials. We developed a Bliss-type analysis to quantify the synergy between anti-PD1 and dabigatran etexilate using published tumor growth data from a mouse model. We then developed a minimal mechanistic computational model to quantitatively study the synergy between anti-PD1 and the thrombin inhibitor dabigatran etexilate in a published study (Metelli et al.) of a preclinical mouse model of colon cancer. Our model included tumor cells, CD8+ T cells, and the pleiotropic cytokine TGF{beta}, whose production in platelets is influenced by thrombin, and described a potential mechanism of interplay among these components in the tumor microenvironment. We performed nonlinear mixed-effects modeling to capture mouse-to-mouse variation in tumor growth under different treatment conditions. The estimated parameter values pointed to several underlying mechanisms of synergy, including increased expansion of CD8+ T cells in the tumor microenvironment in the presence of anti-PD1 and dabigatran etexilate. These predictions can be further validated in future experiments. Thus, the combination of longitudinal tumor growth data, mechanistic population dynamics modeling, and nonlinear mixed-effects modeling can be used to study synergy between ICB and other drugs in controlling tumor growth.

Non-exponential residence times give rise to history-dependent or memory effects in epidemic dynamics, yet their mathematical formulations remain fragmented. Here we develop a unified framework for a history-dependent SEIR model with general latent- and infectious-period distributions. We first establish the equivalence between the stage-age and renewal formulations and show that both admit a unified transition-flow representation. We then introduce the concept of a transition-flow multiplier, which provides a direct link between residence-time distributions and transition flows. When the multiplier is a weighted sum of shifted power functions, the model admits an exact representation in terms of Riemann-Liouville or tempered Riemann-Liouville fractional operators, providing a first-principles foundation for fractional epidemic models. Finally, using Euler-Lotka equations and stochastic-order arguments, we characterize how the latent- and infectious-period distributions shape epidemic growth and decay. A longer latent period slows both early growth and late decay, while, at fixed means, greater variability in the latent period accelerates early growth but slows late decay, whereas greater variability in the infectious period slows both.

Host immunity acquired through prior infections is a key factor shaping the spread of emerging variants. However, accounting for complete infection histories in mathematical models is challenging because the number of possible histories grows exponentially with the number of variants. Here, we present a flexible framework for modeling the co-circulation of multiple variants emerging at arbitrary times while incorporating cross-immunity and variant-specific characteristics. We introduce three tractable approximations: truncation, which retains only infection histories up to a given number of strains; one variant at a time, which is exact when variants circulation does not overlap; and two variants at a time, which is exact when at most two strains co-circulate. We systematically evaluate these approximations across a range of cross-immunity levels and inter-variant emergence delays. Truncation approximations perform best under short emergence delays and strong cross-immunity, whereas the other approaches outperform them at longer delays. When applied to SARS-CoV-2 variant dynamics, the second-order truncation and the two variants at a time approximation accurately capture incidence patterns, whereas the other approaches perform well only for specific subsets of variants. These results identify the parameter regimes in which simplifying assumptions remain valid, provide guidance for model selection, and highlight the key role of infection history in shaping epidemic trajectories.

Devireddy et al. [1] previously showed that an enthalpy-based, macroscale model simulated tissue freezing histories are in close agreement with a coupled model incorporating cellular water transport and intracellular ice formation (IIF). That comparison, however, was performed for a limited set of tissue geometries, cooling conditions, and biophysical parameters. Here, we revisit that conclusion by reimplementing the published model from its governing equations and systematically expanding the parameter space to determine when microscale processes alter the predicted macroscale freezing response. We incorporate both surface catalyzed nucleation (SCN) and volume catalyzed nucleation (VCN) of intracellular ice with cellular osmotic dehydration and examine a broad range of tissue dimensions, cooling rates, convective boundary conditions, membrane permeabilities, cell sizes, activation energies, intracellular water fractions, and nucleation rates. Linearization of the cellular water transport equation about osmotic equilibrium yields an osmotic relaxation timethat can be compared with the local residence time of tissue undergoing phase change. Their ratio defines a dimensionless group, that organizes the transition between regimes in which the uncoupled enthalpy formulation and the coupled micro and macroscale formulation produce similar or substantially different results. Across the one-dimensional parameter space examined here, the two models remain within a prescribed thermal-error bound for the dimensionless group below approximately 0.1, whereas the discrepancy increases systematically above this threshold. Model divergence also depends on intracellular ice nucleation, i.e., when nucleation is sufficiently slow, cellular water transport can become rate limiting for latent heat release. Conversely, sufficiently rapid intracellular ice formation reduces the macroscopic discrepancy by providing an additional pathway for latent-heat release. Multidimensional simulations further show that geometry alone can move a nominally conventional freezing protocol into the coupled regime. In particular, multidirectional cooling can reduce the local residence time sufficiently to produce substantial differences in predicted frozen volume even when the corresponding thermal histories appear similar. Finally, we formulate an inverse approach for estimating IIF parameters from tissue scale thermal histories. Synthetic simulations indicate that parameter recovery is feasible within a finite range of cooling conditions, providing a testable strategy for future experimental measurements. Together, these results provide a mechanistic criterion for deciding when microscale biophysical processes must be retained in macroscale tissue-freezing models.